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(57) ABSTRACT 

Real time battery impedance spectrum is acquired using one 
time record, Compensated Synchronous Detection (CSD). 
This parallel method enables battery diagnostics. The excita- 
tion current to a test battery is a sum of equal amplitude sin 
waves of a few frequencies spread over range of interest. The 
time profile of this signal has duration that is a few periods of 
the lowest frequency. The voltage response of the battery, 
average deleted, is the impedance of the battery in the time 
domain. Since the excitation frequencies are known, synchro- 
nous detection processes the time record and each compo- 
nent, both magnitude and phase, is obtained. For compensa- 
tion, the components, except the one of interest, are 
reassembled in the time domain. The resulting signal is sub- 
tracted from the original signal and the component of interest 
is synchronously detected. This process is repeated for each 
component. 
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METHOD OF DETECTING SYSTEM 
FUNCTION BY MEASURING FREQUENCY 
RESPONSE 

CROSS-REFERENCE TO RELATED 5 

APPLICATIONS 

This application is a continuation of U.S. patent applica- 
tion Ser. No. 11/313,546, filed Dec. 20, 2005, now aban- 
doned, which claims the benefits of U.S. Provisional Patent to 
Application Nos. 60/637,969, filed Dec. 20, 2004 and 60/724, 
631, filed Oct. 7, 2005. The disclosures of each of these 
applications are hereby incorporated by reference in their 
entirety, including all figures, tables and drawings. 

The subject invention was made with government support 15 
under a research project supported by NASA, Grant No. 
NNA05AC24C. The government has certain rights in this 
invention. 

BACKGROUND OF THE INVENTION 20 

Electrochemical Impedance Measurement Systems use the 
Bode analysis technique to characterize the impedance of an 
electrochemical process. It is a well established and proven 
technique. The battery being evaluated is excited with a cur- 25 
rent that is single frequency and its response is measured. The 
process is repeated over a range of frequencies of interest 
until the spectrum of the impedance is obtained. The method 
is effective but time consuming, as the process is serial. A 
parallel approach using band width limited noise as an exci- 30 
tation current can obtain the same information in less time. 
The system response to the noise is processed via correlation 
and Fast Fourier Transform (FFT) algorithms and many such 
responses are averaged. The result is the spectrum of response 
over the desired frequency range. The averaging of many 35 
responses also makes this process somewhat serial. Another 
technique assembles the current noise waveform from a sum 
of sinusoids each at a different frequency. The system 
response as a time record is acquired and processed with the 
FFT algorithm. To reduce noise multiple time records of 40 
waveforms are processed and their resultant spectra averaged. 
This process is also serial. 

There remains a need for real time acquisition of battery 
impedance for control and diagnostics over a limited fre- 
quency range. This method of acquisition should be a true 45 
parallel approach that uses a single time record of battery 
response with a duration compatible to a real time control 
process. 

SUMMARY OF THE INVENTION 50 

The invention involves using a parallel approach to analyze 
battery impedance or other system functions. A number of 
frequencies are selected over which the battery is to be tested. 
These frequencies are assembled into an excitation time 55 
record (ETR) that is the sum of the sinusoids of the frequen- 
cies and the length of such periods of the lowest of the fre- 
quencies. The ETR is conditioned to be compatible with the 
battery. The battery is then excited with the ETR and the 
response time record (RTR) is captured. The RTR is then 60 
synchronized to the ETR and processed by a series of equa- 
tions to obtain a corrected time record (CTR). The CTR is 
reassembled to obtain an estimated corrected time record 
(ECTR). The ECTR is subtracted from the CTR to get an 
error. The error is minimized to achieve the frequency 65 
response estimate. Error is minimized using compensated 
synchronous detection (CSD) using a CSD algorithm which 
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can be implemented by a neural network. The subject method 
allows a parallel implementation for swept frequency mea- 
surements to be made utilizing a composite signal of a single 
time record that greatly reduces testing time without a sig- 
nificant loss of accuracy. 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 shows the Filter Ideal Magnitude Response. 

FIG. 2 shows the Filter Ideal Phase Response. 

FIG. 3 shows the Filter Uncompensated Synchronous 
Detected Magnitude Response. 

FIG. 4 shows the Filtered Compensated Synchronous 
Detection (CSD) Magnitude Response. 

FIG. 5 shows the Filter CSD Phase Response. 

FIG. 6 shows the Lumped Parameter Model (LPM). 

FIG. 7 shows a portion of the Sum of Sines (SOS) signal to 
LPM 13 lines 10 periods. 

FIG. 8 shows a portion of the LPM time response 13 lines 
10 periods. 

FIG. 9 shows the LPM ideal magnitude response 13 lines 
10 periods. 

FIG. 10 shows the LPM ideal phase response 13 lines 10 
periods. 

FIG. 11 shows the LPM uncompensated magnitude 
response 13 lines 10 periods. 

FIG. 12 shows the LPM CSD magnitude response 13 lines 
10 periods. 

FIG. 13 shows the LPM CSD phase response 13 lines 10 
periods. 

FIG. 14 shows a portion of the LPM SOS signal 25 lines 10 
periods. 

FIG. 15 shows a portion of the LPM time response 25 lines 
10 periods. 

FIG. 16 shows the LPM ideal magnitude response 25 lines 
10 periods. 

FIG. 17 shows the LPM ideal phase response 25 lines 10 
periods. 

FIG. 18 shows the LPM uncompensated magnitude 
response 25 lines 10 periods. 

FIG. 19 shows the LPM CSD magnitude response 25 lines 
10 periods. 

FIG. 20 shows the LPM CSD phase response 25 lines 10 
periods. 

FIG. 21 shows the Mean Squared Error (MSE) comparison 
for low pass filter frequency response. 

FIG. 22 shows the MSE comparison for detection of LPM 
impedance. 

FIG. 23 is a flow chart showing the method of the subject 
invention. 

DETAILED DESCRIPTION OF THE INVENTION 

A system function can be identified over a limited number 
of specific frequencies. The desired frequencies are 
assembled as an excitation time record that is a sum of those 
sinusoids and have a length of several periods of the lowest 
frequency. The time step selected must be compatible with 
Shannon’s sampling constraints for the highest frequency 
component. The individual waveforms should be sine waves 
of equal amplitude but alternating signs for a phase shift of 
1 80 degrees between the components. Alternating 1 80 degree 
phase shift will minimize a start up transient. The RMS and 
the rogue wave peak (sum of the absolute values of all com- 
ponent peaks) of the assembled time record must be compat- 
ible with the system being excited and the Data Acquisition 
System (DAS) that will capture the response. 
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The Excitation Time Record (ETR), as a current for imped- 
ance identification or voltage for system function identifica- 
tion, is signal conditioned to be compatible with the Unit 
Under Test (UUT). As part of the signal conditioning, anti- 
aliasing filters insure that all frequencies generated by the 5 
digital to analog conversion process other than the intended 
frequencies are suppressed. The UUT is excited by the ETR 
and a time record of the UUT response is captured by the 
DAS. The UUT response time record (URTR) is synchro- 
nized to and the same length as the ETR. to 

In order to fit steady state sinusoidal response assumptions 
a pre-selected number of data points R at the beginning of the 
URTR must be discarded. In general the sum of those data 
points total to a time that is larger than the transient response 
time of the UUT at the front end of the ETR. The UUT 15 
corrected time response is referred to as the CTR. 

The first estimate of components, magnitude and phase, of 
the frequency response is made by processing the URTR via 
Equations 3 through 10. It is important that a zero mean be 
established prior to processing with Equations 3 to 10 (see 20 
below). 

The core of this whole concept is that an estimate of the 
UUT corrected time response, the CRT is made by reassem- 
bling the URTR using the estimates of the individual fre- 25 
quency components with the same time step and then discard- 
ing the first R time steps to become the ECRT. The difference 
between the CRT and the ECRT is an error and minimizing 
this error will increase the accuracy of the frequency response 
estimates. 30 

A first approach to minimizing the error between the CRT 
and the ECRT is Compensated Synchronous Detection 
(CSD). The CSD algorithm synthesizes a residual time record 
of the original time record using the magnitudes of the in 
phase and quadrature components for each frequency except 35 
the one to be detected. This synthesized residual is then sub- 
tracted from the original time record. The resulting compen- 
sated time record is processed with synchronous detection 
and a new compensated estimate of the response at the detec- 
tion frequency is obtained. Since all of the other components 40 
in this compensated time record are suppressed the error from 
leakage at those other frequencies is less. This process is 
repeated for each of the frequencies. Assembling the residual 
time record and generating the compensated time record are 
illustrated by Equations 1 1 and 12 (see below). 45 

Another approach to minimize the error between the CRT 
and the ECRT is to use a Neural Network. The first estimate 
of the component magnitudes and phases is made as for the 
CSD technique. Those values are stored and the ECRT is 
calculated. This signal is then subtracted from the CRT to 50 
produce a response residual. The synchronous detection is 
then performed upon this residual and the component mag- 
nitudes and phases are again stored. These components are 
then used to reconstruct an estimate of the residual signal. 
This estimated residual signal is subtracted from the initial 55 
residual signal to produce a residual, residual signal. This is 
then synchronously detected and the loop starts again. This is 
repeated as many times as desired, each time with the result- 
ant components being stored. The assumption is that there is 
a functional relationship between these resultant components 60 
and the true system response. A neural network is then used to 
determine this relationship. The previously stored results 
become the test dataset for the neural network. The network 
has been trained previously on a similar unit and a known 
ideal response (e.g. battery impedance measured using elec- 65 
trochemical impedance spectroscopy). The output of the net- 
work is the estimate of the response. 
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A further step of the subject method is to shift the complete 
set of desired frequencies and repeat the whole process. This 
step could be repeated many times with different shifts to 
develop a high resolution frequency response that for battery 
impedance could be comparable to that provided by electro- 
chemical impedance spectroscopy. Thus the subject method 
can provide capability for both limited frequency response in 
real time or high resolution frequency response not in real 
time for periodic system in depth diagnostics. 

The system of the subject invention is based on the follow- 
ing theoretical design. The unit under test is excited with a 
limited sum of sinusoids each at a different frequency that is 
spread over the range of interest. The magnitude, frequency 
and phase of each sinusoid making up the sum are known. If 
the total response of the system is measured via a sample data 
system at an acceptable sample rate and an adequate duration 
time record is acquired, then a simple algorithm that uses the 
known magnitude, frequency and phase of each individual 
sinusoid will process the single time record. This analysis will 
obtain the true Bode response at the selected frequencies 
spread over the range of interest all in parallel. The following 
synchronous detection analysis is the basis of this simple 
algorithm. The reference waveform is chosen as a sine as at 
time zero everything will be at zero. 

Equation 1 gives relationship for the parallel excitation. 


Ut) = ^ A,sin(fc>,r + 0in, ) 


Equation 2 gives the measured sampled data response of 
the system 


foutU] = ^ B}Sm((x)i(j- 1)A t + <pouti)\ j = 1 :N 


Where: 

A z amplitude of the \ th input sinusoid 
B z amplitude response of the \ th output sinusoid 
co z radian frequency of the \ th sinusoid 
At time step of data system 
(j)in z phase of the 1 th input sinusoid 
(|)out z phase response of the \ th output sinusoid 
N number of points of the response time record 
M number of different sinusoids of the excitation time 
record 

Each component magnitude and phase of the system 
response at all the excitation frequencies can be obtained via 
the following synchronous detection analysis. This analysis 
quantifies the response at the k th radian frequency c% with the 
“in phase” and “quadrature” response. The analysis incorpo- 
rates the feature of discarding a pre-selected number of points 
R at the beginning of the system response in order to meet the 
assumption of steady state sinusoidal response. Additionally, 
for most applications prior to processing the data the mean of 
the acquired time record should be computed and deleted. 
The presence of a non-zero mean could corrupt the estimate 
of the lowest frequency component. 
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Equations 9 and 10 give the final results for synchronous 
detection. This process is repeated for all M of the excitation 
sinusoids. As stated, the results depend on the time record 

(3) 

being infinite in duration. This is not the case and thus the 
5 results are in error by leakage from the other components. 
This error can be reduced without increasing the time record 
using the following algorithm for compensated synchronous 
detection (CDS). 

The CSD algorithm synthesizes a residual time record of 
the original time record using the magnitudes of the in phase 
and quadrature components for each frequency except the one 
to be detected. This synthesized residual is then subtracted 
from the original time record. The resulting compensated 
time record is then processed with synchronous detection and 
a new compensated estimate of the response at the detection 
frequency is obtained. Since all of the other components in 
this compensated time record are suppressed the error from 
leakage at those other frequencies will be less. This process is 
repeated for each of the frequencies. Assembling the residual 
(5) 20 time record and generating the compensated time record are 
illustrated by Equations 11 and 12. 
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15 


If the time record were infinite the summation over j would 
average everything to zero except the final result of Equation 25 
5. The quadrature analysis follows in the same 

Quadrature: 


IrkUI- Yj (j- 

p=\,p*k 


l)Ar) + F Qp cos(co p (j - l)Ar)); 
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j=R + 1 


AKCos(a)K(j - l)Ar + 0inft)^ 5,sin(w,(jf - l)Ar + pouti ) 
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Where: 

^Out is the original time record 
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( 8 ) 


Again the summation over j for infinite time record aver- 50 
ages everything to zero except the final result of Equation 8. 
Equations 5 and 8 can be combined to give magnitude and 
phase for the k th frequency response. 


\Foutfr | — yj foutf' fy‘t 


, 2 

'OUtfc 


(9) 


AkBk I I ~ A^Bk 

2 y sin z (pm K - pout K ) + cos 2 (pm K - pout K ) = — — 


CfK 0w , is the compensated time record 
F p estimated in phase amplitude response at the p tf 


55 


f fre- 
quency 

F Qp estimated quadrature amplitude response at the 
frequency 

co > radian frequency of the p th sinusoid 
At time step of data system 
N number of points of the output time record 
M number different sinusoids of the excitation function 
R number of points of the output time record that are 
discarded 

The synchronous detection algorithm described by Equa- 
tions 1 through 8 is applied to the compensated time record of 
Equation 12 and better estimates of and F ^ are obtained. 

This process can be repeated again until the total difference 
between a completely synthesized time record response and 
the original time record is minimized. 

The following examples are offered to merely illustrate the 
method of the subject invention and should not be construed 
as limiting. 

EXAMPLE 1 

Analytical Testing on a Sum of Sines 


□F™ 




sin(d>outv — (b iiuO 

(10) 60 

( Fq out K j 

= tan 1 


= 

l Fout K J 

AkBk 

K COS (poutK - P m K.) ^ 


( sin (poutR ~ 

tan — 

Vcos {poutK ■ 


pin K ) > 


0m/OJ 


= (pout K ~ 0 in*) 


65 


The CSD algorithm was evaluated using a simple signal 
that was assembled from a finite sum of equal amplitude sin 
waves (Sum of Sines, SOS) with frequencies distributed loga- 
rithmically over a limited range. The objective of the analysis 
was to assess how well the CSD algorithm could pick out the 
amplitude for each component. 

To check out the concept analytically a MATLAB matrix 
calculation computer software code was written that was a 
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logarithmic mix of 5 equal unity amplitude frequencies (5 0 ' 5 , 
5 1 , 5 1 ' 5 , 5 2 , 5 2 5 Hz). The acquired time record was set to 10 
periods of the lowest frequency and the time step was set to 
Yio of the period of the highest frequency. As per Equation 9, 
error- free detection should estimate the amplitude of 0.5 for 
each component. Table 1 gives the estimate for the first pass or 
simple synchronous detection and the second pass, the CSD 
algorithm. The MATLAB matrix calculation computer soft- 
ware code for this analysis is given in Morrison et al., 2005. 


TABLE 1 


Compensation algorithm analytical results 


1st Pass (Simple Synchronous Detection) 


Frequency 

5 5 

5 1 

5 i . 5 

5 2 

5 25 

Amplitude 

0.5060 

0.5060 

0.4975 

0.4988 

0.5008 


2nd Pass (CSD Algorithm) 


Frequency 

5 5 

5 1 

5 1-5 

5 2 

5 2 ' 5 

Amplitude 

0.5004 

0.4998 

0.5004 

0.5000 

0.5003 


As seen in Table 1 analytically, the compensation tech- 
nique does appear to work as the error for the second pass is 
much reduced. This initial check of the concept was applied to 
detect amplitude only and no phase detection. The signal is a 
sum of equal amplitude sin waves being decomposed into the 
individual components by the algorithm. 

EXAMPLE 2 

Analytical Testing of a Low Pass Filter 

A recursive model of a second order low pass function was 
excited with a SOS input signal. The CSD algorithm was then 
used to estimate the frequency response at each of the specific 
frequencies making up the SOS. 

A spread of 13 specific frequencies was chosen that were 
spaced in l A decade steps starting from 0.1 Hz up to 100 Hz. 
Using these frequencies a mix of equal unity amplitude sine 
waves was created. This range of frequencies was picked as 
the research performed at the Idaho National Laboratory with 
batteries is over this frequency spread. The signal was dis- 
cretized with a time step that was 10% of the period of the 
highest frequency. The length of the time record was set at 1 0 
periods of the lowest frequency. 

A recursive model of a second order Butterworth low pass 
function was developed. The center frequency was set at the 
middle of the SOS frequency spread. The filter response to the 
SOS input time profile was computed. 

Finally, the CSD code was developed to estimate the filter 
frequency response from the time profile generated by the 
recursive model code. That code is the implementation of 
Equations 1 through 12. Additionally, the implementation has 
the ability to discard a number of user selected points at the 
beginning of the time profile such that the remaining data 
better fits the assumption of “steady state” response. The 
following 6 figures are MATLAB matrix calculation com- 
puter software discrete plots of the ideal frequency response, 
the uncompensated frequency response and finally, the com- 
pensated frequency response. All 3 have magnitude plotted 
against the Log of frequency. The FIGS. 1 and 2 are the ideal 
magnitude and phase response. FIG. 3 is the uncompensated 
magnitude response FIG. 4 is the compensated magnitude 
response. The improvement seen by comparing the last two 
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figures is obvious. FIG. 5 is the compensated phase response. 
Error at the higher frequencies for the compensated result is 
likely due to processing a small signal at the high frequencies 
relative to larger signals at the lower frequencies. The net 
5 result i s that one time record yields a limited number of points 
for both magnitude and phase. This technique shows promise 
for real time applications. The MATLAB matrix calculation 
computer software code for this analysis is given in Morrison 
et al., 2005. The next approach will apply the concept ana- 
l o lytically in an attempt to identity the impedance of a computer 
model of a battery. 

EXAMPLE 3 

15 Analytical Testing of CSD with a Battery Model 

The CSD algorithm was evaluated analytically via a com- 
puter simulation of the detection of the impedance of the 
Lumped Parameter Model of a battery (LPM) that was devel- 
20 oped by the Idaho National Laboratory (INL) (see, Freedom- 
Car Battery Test Manual, 2003). A computer model for the 
LPM that will simulate battery voltage response to an arbi- 
trary battery current was also developed at INL by Fenton et 
al. 2005 . The voltage response of the model normalized to the 
25 current in the frequency domain will be the battery imped- 
ance. The equivalent circuit for the LPM with parameter 
identification is given by FIG. 6. 

The LPM was excited with a current source \ IN that was a 
SOS and the CSD algorithm was used to identity the imped- 
30 ance seen looking into the LPM over a limited range of 
discrete frequencies. It should be noted that the polarity of the 
voltage response was defined as negative because the SOS 
excitation current was a discharge (negative relative to imped- 
ance). The CSD algorithm was used to obtain the frequency 
35 response of the LPM. The CSD response magnitude and 
phase were compared to the ideal response. Initially, the 
algorithm failed to match ideal impedance. The response of 
the battery terminal voltage to a discharge SOS current signal 
will contain a DC term caused by the DC battery voltage. 
40 Synchronous detection of any specific frequency in the 
response will cause a noise frequency in the resultant spec- 
trum at the frequency being detected and, at amplitude of the 
product of the DC battery voltage and the amplitude of the 
detection signal. That noise amplitude will be large relative to 
45 the signal being detected. Averaging enough cycles in the 
resultant time record will reject this noise. However, for a real 
time application the length of the time record needs to be short 
and not long. The problem was resolved when the mean was 
deleted from the prediction response of the LPM. The number 
50 of frequency lines was set to 1 3 and they were logarithmically 
spread from 0.01 Hz to 1 .0 Hz. The time record was set to 10 
periods of the lowest frequency. In the CSD algorithm no data 
points were discarded. Table 2 gives the analysis specifics 
with LPM data that is typical for some Li batteries that INL 
55 had recently tested. INL performed the testing per methods in 
FreedomCar Battery Test Manuals 2003. 

TABLE 2 


6Q Representative LPM and Analysis Data 

Voc = 3.8 

Cp = 666.6667 At Rest 

Coc = 1.6667e+003 At Rest 

Ro = 0.0250 
Rp = .0150 

65 M=13 Number of frequency lines 

Dt = .01 Time step, sec 
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TABLE 2 -continued 



Representative LPM and Analysis Data 

N = 100000 

Total number of points 

F = .01 

Starting Frequency, Hz 

FF = 10 

Frequency spread in decades 

S = .125 

Step size (log) over the decades, 8 steps per 
decade 

NN= 10 

Length of time record in number of periods 
of lowest frequency 


The following 7 figures are the plots of the analysis results. 
FIG. 7 is the time record of the SOS current signal. FIG. 8 is 
the time record of the LPM voltage response. FIGS. 9 and 10 
are the LPM ideal impedance magnitude and phase. FIGS. 11 
and 12 are the uncompensated and the compensated magni- 
tude response. FIG. 13 is the compensated phase response. 

The number of frequency lines was increased to 25 and 
everything else left the same. The following 7 figures illus- 
trate the results. FIG. 14 is the time record of the SOS current 
signal. FIG. 15 is the time record of the LPM voltage 
response. FIGS. 16 and 17 are the LPM ideal impedance 
magnitude and phase. FIGS. 18 and 19 are the uncompen- 
sated and the compensated magnitude response. FIG. 20 is 
the compensated phase response. 

Additional cases run showed that as the length of the 
acquired time record in the number of periods of the lowest 
frequency gets cut back, the number of frequency lines that 
can be resolved without a big increase in error must be cut 
back. For example, 5 periods with 25 lines gave terrible 
results but 5 periods with 13 lines was fine. These positive 
results are only analytical. Nevertheless, they offer promise of 
positive expected performance when applied to a physical 
system. All the MATLAB matrix calculation computer soft- 
ware codes for this analysis are given by Morrison 2005. 

EXAMPLE 5 

Neural Network Enhanced Synchronous Detection 

In order to improve the accuracy of the CSD, studies were 
conducted upon neural network enhancement of the detection 
of the individual frequency components of the response of a 
linear system to a SOS input signal. The concept is very 
similar to the CSD, with some slight changes and a neural 
network output layer. For a second order lowpass filter, ordi- 
nary synchronous detection performed on the filter response 
showed a mean squared error (MSE) of 2.6x1 0 -3 . The com- 
pensated synchronous detection technique displayed a MSE 
of 1.6xl0 -3 . The neural network enhanced synchronous 
detection showed a MSE of 0.2xl0 -3 . Results for the lumped 
parameter model of a lithium ion battery were similar. 

The theory of the neural network enhanced synchronous 
detection is based on the classical synchronous detection and 
the inherent error associated with time records of finite 
length. 

Given an input signal comprised of a sum of N sinusoids: 


N 

x(t) = ^ x{OJ n t) 

n = 1 


10 

the output of the system would be: 


I V 

5 K0 = Yi y ^ nt ^ 

n= 1 

Being that the input and output are sinusoids, they are 
assumed to have started at t=-oo and continues to t=oo. In 
10 reality this is not the case, however. The time record of the 
signal is first of all finite in length and second of all, it is 
sampled. Given a sampling frequency of at least the Nyquist 
frequency (twice the highest frequency in the signal), the 
signal can be reconstructed without error. This does not, 
15 however, rectify the finite time length of the signal. Because 
of this, errors enter into the synchronously detected frequency 
components. 

Synchronous detection involves multiplying the acquired 
signal by sines and cosines of the desired frequencies and 
20 summing the results, as shown below: 


1 

®oj n = lim^oo Jilu 
j=~i 

1 V 1 

fiion = 2iZj yificosRw'Ar] 


3° 

The magnitude of a given frequency is obtained as follows: 

M UJft =2yja (0ft 2 +^ UJft 2 

and the phase may also be obtained as follows: 

35 



40 

where: 

y[j]=the sampled signal, 

a=the sine component of the response for the frequency a>„ 
45 p=the cosine component of the response for the frequency 

av 

At=the sample time step. 

Notice that the summations are infinite in length. In appli- 
cation, the summation would be from 0 to the length of the 
50 recorded signal . If the time record was infinite, then any errors 
would cancel out. Since our time record is not infinite, errors 
remain in the calculated response. 

It is important to note that an estimate of the original signal 
can be reconstructed by summing the sine and cosine signals 
55 of the different frequencies where each sine signal is multi- 
plied by its sine response component and each cosine signal is 
multiplied by its cosine response component. The resulting 
estimate would be exact if the time record was infinite, but 
since it is of finite length, the reconstructed signal isn’t 
60 exactly correct. This may be viewed as containing noise. 


N 

y[i 1 = Y ^sin(w„/Ar) + A, n cos(w„/Ar) + rj [z] 

n=l 


65 



11 


US 7,395,163 B1 


where: 

r|(i)=the noise due to error 

At=the sample time step 

i=the ith position in the acquired time record 5 

N=the number of frequency lines 

In the compensated synchronous detection approach, this 
difference between the actual signal and the reconstructed 
signal is exploited to filter out all but one of the frequencies 
and increase the accuracy of the measurement (see above) .In 10 
Neural Network Enhanced Synchronous Detection 
(NNESD), the approach is slightly different. The premise is 
that the residual left from subtracting the reconstructed signal 
from the actual signal still contains some useful information. 

r[i]=y[i]-y[i] 

For each frequency of interest, the sine and cosine compo- 
nents are detected on the original signal. These are then used 
to reconstruct the signal and then stored. The reconstructed 20 
signal is then subtracted from the original signal to obtain the 
residual signal. Synchronous detection is then performed 
upon this residual signal; once again the sine and cosine 
components are used to reconstruct the residual signal and 
also stored. This second residual signal is then subtracted 25 
from the first residual signal. This loop continues until a 
sufficient number (M) of sine and cosine responses are 
obtained for each frequency. Now the assumption that is made 
is that there is some functional relationship between these M 
responses and the true responses. 30 

a oj„2? (U„2> • ■ • > a uyi/5 a ayi/) 

IU„1? a uj„2? Poj„ 2> • ■ ■ 5 a aj„M> a tD^f) 

This functional relationship most likely differs from sys- 35 
tern to system and also based upon system operating condi- 
tions. To deduce what this relationship is for any system 
would be time consuming. Instead, a generalized regression 
neural network is implemented to predict the true response. 

A Generalized Regression Neural Network (GRNN) is a 40 
form of a radial basis neural network suitable for problems 
involving regression, function estimation, or prediction of a 
continuous value Wasserman, 1993, Alpaydin, 2004. This is 
in contrast to a very similar network, the Probabilistic Neural 
Network, which is used for classification applications. Unlike 45 
multi-layer perceptrons, which require training, a GRNN is 
constructed from a training set of example inputs and the 
corresponding outputs. The spread of the radial basis function 
is used to compute the bias. The spread in effect controls the 
spread of the radial basis function. The example inputs are in 50 
the form of an mxn matrix where each of the m rows repre- 
sents an observation and each of the n columns represents a 
feature, or parameter. The corresponding example output is 
an mxl vector. The input is said to be n-dimensional. The 
network is divided into 2 layers. The first, or input layer, 55 
consists of the example inputs. The output layer consists of 
the example outputs. The network generates its output in the 
following manner. The geometric distance is calculated 
between the newly presented input and each of the example 
inputs in the input layer: 60 


d m 



(x n -IW n ) 2 


65 
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This produces a vector of length m. Each element of the 
vector is then multiplied by the bias, which is 0.8326/Spread. 
The vector then goes through the radial basis function. 


The radial basis function produces an output that gets closer 
to 1 as the input gets closer to 0. The resulting vector is a 
ranking of how close each of the example cases are to the new 
input. The normalized dot product is then calculated between 
the vector and the example output vector. This is the network 
output. 


v = — V x m LW m 
y m m 

m=l 


The NNESD is based upon the GRNN approach. Analyti- 
cal testing will provide preliminary validation of the concept. 

The NNESD concept was analytically tested on a 2 nd order 
Butterworth filter using the MATLAB matrix calculation 
computer software BUTTER function and the INL LPM 
FreedomCar Battery Test Manual, 2003, Fenton et al., 2005 
for the lithium ion battery. An SOS input signal was used for 
both cases. For the filter, the component frequencies were 
varied for each run in order to build up a training set with 10 
component lines to be detected. The filter was run 1 00 times 
and the output and target response was calculated. The data 
was randomized and half of the data was used to construct the 
GRNN and the other half was used to verify the network. The 
mean squared error (MSE) of the predicted value was calcu- 
lated and the process was repeated 20 times, shuffling the data 
each time. The results are shown in FIG. 21. The MSE of both 
the standard synchronous detection and the CSD are also 
shown for comparison. The MATLAB matrix calculation 
computer software code used for this analysis is given in 
Morrison, 2005. 

The technique was then tested on the LPM FreedomCar 
Battery Test Manual, 2003, Fenton et al., 2005. The input 
parameters to the model were nominally set as per Table 2 and 
randomly varied by up to 5% each run for 1 00 runs to generate 
the dataset. The number of lines, frequency spread and time 
step was held as per Table 2. The same training and testing 
scheme that was outlined above was used. The results are 
shown in FIG. 22. The MATLAB matrix calculation com- 
puter software code used for this analysis is given in Morri- 
son, 2005. 

This limited analytical validation has shown that the 
NNESD concept will significantly reduce the error in the 
estimate of the frequency components of a given system 
response signal. This concept allows a parallel implementa- 
tion for swept frequency measurements to be made utilizing a 
composite signal of a single time record that greatly reduces 
testing time without a significant loss of accuracy. 

CONCLUSIONS 

The physical validation of the CSD or NNESD concept 
will rely heavily on the work performed by W. Albrecht in his 
thesis research: “Battery Complex Impedance Identification 
with Random Signal Techniques” Albrecht, 2005. In this 
approach, a National Instruments data acquisition and pro- 
cessing system was used along with a custom analog condi- 
tioning system. The CSD or NNESD algorithm could be 
installed directly on that system. The system software will be 
CSD or NNESD rather than Noise Identification of Battery 
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Impedance (NIBI). The NIBI approach would acquire about 
100 time records of the battery response to noise current. 
Clearly, a time record would have to be of a length of multiple 
periods of the lowest frequency of interest. The CSD or 
NNESD approach will acquire one time record of a length of 5 
multiple periods (exact number is still to-be-determined) of 
the lowest frequency of interest. The excitation signal would 
be generated analytically with software as in the NIBI system. 
The analytical signal would be preconditioned with a digital 
low pass filter as in the NIBI system. The CSD or NNESD to 
system may require an analog filter prior to the current driver 
and after the D/A. This analog filter at the current driver could 
serve as the prime anti-aliasing filter. This system will use the 
same bias compensation approach to remove most of the DC 
battery voltage from the acquired signal. Improved noise 15 
rejection and increased sensitivity could be achieved if the 
voltage sensing were upgrade to full differential via a 4-wire 
system rather than the 2-wire single-ended system the NIBI 
used. The increase in sensitivity will enable a reduction in the 
level of the excitation signal required. It is anticipated that the 20 
sampled voltage will be processed directly with the CSD or 
NNESD algorithm. It is also anticipated that a system cali- 
bration would be done exactly as the NIBI system by mea- 
suring the impedance of the test leads to determine any sys- 
tem measurement offset and phase shift. 25 

It is understood that the foregoing examples are merely 
illustrative of the present invention. Certain modifications of 
the articles and/or methods employed may be made and still 
achieve the objectives of the invention. Such modifications 
are contemplated as within the scope of the claimed inven- 30 
tion. 
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The invention claimed is: 

1. A method for detecting system function of a unit under 
test by measuring frequency response, the method compris- 60 
ing the steps of: 

(a) selecting a number of frequencies over which the func- 
tion of the unit under test will be tested; 

(b) assembling an excitation time record that is the sum of 
the sinusoids of the frequencies and a duration of greater 65 
than or equal to a period of the lowest selected fre- 
quency; 


(c) conditioning the excitation time record to be compat- 
ible with the unit under test; 

(d) exciting the unit under test with the excitation time 
record and simultaneously capturing a response time 
record with a data acquisition system; 

(e) processing the response time record using the following 
synchronous detection equations to obtain estimated 
frequency components of magnitude and phase for one 
of the selected frequencies: 


" N-R 

N 


^ - l)Ar + 0injr)^ Bism(u>i(j - l)Ar + <pouti ) 


c 1 V (AkBk, 
°“' K ~ N-R Z . J I 2 1 

j=R+\ ' 


cos(0injf — (poutK ) — cos(2oj k{J - l)Ar + <pm^ + 0owr^)] h 
A K B ir 

i±k = 1 


A k B{ 

2j ~ 2 ~ [cos((w* - QJi)(j - l)Ar + <pm K - <pouti) - 


cos((w* + ojj){j - l)Ar + <pmK + (pouti)] 


A K B K 

F ou t K - — — cos(0inK -<pout K ) 


F %utr = 


1 


' ut K N - R 


\ < Af(Cos(o)fc(j - l)Ar + 0in^)^ - l)Ar + <pouti) > 


F 1 Y \A K B K 

q oiut K N _ R 1 2 L 

j=R + 1 ^ 

sin (<poutx — + sin(2 (Ox(j — l)Af + 0in^ + <poutx)] + 

Vi A K Bi 

2j ~ 2 ~ [sin((^A: + o>i)(j ~ l)Af + <pm K + <pout t ) - 


sin((<t>* - - l)Ar + 0in^ - <pouti)] 
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Fq outK = sin {<pout K - <pm K ) 
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Vcos {(poutK - ) 


Where: 

A z amplitude of the \ th input sinusoid 

B z . amplitude response of the \ th output sinusoid 

oo z radian frequency of the \ th sinusoid 

At time step of data system 

4>in z phase of the \ th input sinusoid 

(|)out z phase response of the \ th output sinusoid 
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N number of points of the response time record Where: 

M number different sinusoids of the excitation time record f is the original time record 

R the number of points in the response time record that may f j s correction time record 


be discarded 

K is the frequency index for the system function being 
detected; 

(f) repeating step (e) to obtain estimated frequency com- 
ponents for each selected frequency; 

(g) assembling the estimated frequency components to get 
an estimated response time record; 

(h) subtracting the estimated response time record from the 
captured response time record to get an error; and 

(i) minimizing the error to achieve the frequency response. 

2. The method of claim 1, wherein said error is minimized 

using the following equations which create a corrected time 
record and a compensated time record: 


Irk [j] = Yj U - Ar ) + f Qp C 0 ^pU ~ 1 ) A? )); 

P=l,P*K 

j = R+l:N 


QKoJjMoJjHrkD] 


Cf KOut is the compensated time record 
F estimated in phase amplitude response at the p th fre- 
quency 

F Qp estimated quadrature amplitude response at the p th 
frequency 

io co radian frequency of the p th sinusoid 
At time step of data system 
N number of points of the response time record 
M number of different sinusoids of the excitation time 
record 

1 5 R number of points of the response time record that may be 

discarded 

K is the frequency index for the system function being 
detected. 

20 3 . The method of claim 1, wherein said error is minimized 

using a neural network. 

4. The method of claim 1, wherein said unit under test is a 
battery. 



